Differential Equations Solving
Michael Liou
4/20/2022
Software
In R, there are a number of packages that deal with differential equations. CRAN View
deSolve- the main integratorReacTran- simulating reaction diffusion equationsFMA- “flexible modelingbvpSolve- boundary value problemssimecol- interactive enviromnent for implementing models
Stochastic Differential Equations have a different set
sde- simulatio n and inference of stochastic differential equationsGillespieSSAadaptivetau- are for quick simulation
deSolve Package
FME Package
Stands for flexible modeling pacakge, allows one to define some cost function
rootsolver
ODE
Using ode from package deSolve
Example: Prey and Predator
LVmod <- function(Time, State, Pars) {
with(as.list(c(State, Pars)), {
Ingestion <- rIng * Prey * Predator # \beta x y
GrowthPrey <- rGrow * Prey * (1 - Prey/K) # \alpha x
MortPredator <- rMort * Predator
dPrey <- GrowthPrey - Ingestion
dPredator <- Ingestion * assEff - MortPredator
return(list(c(dPrey, dPredator)))
})
}
pars <- c(rIng = 0.2, # /day, rate of ingestion
rGrow = 1.0, # /day, growth rate of prey
rMort = 0.2 , # /day, mortality rate of predator
assEff = 0.5, # -, assimilation efficiency
K = 10) # mmol/m3, carrying capacity
yini <- c(Prey = 1, Predator = 2) # initial states of the system
times <- seq(0, 200, by = 1) # initial time
out <- ode(yini, times, LVmod, pars)
summary(out)## Prey Predator
## Min. 1.0000000 1.8632829
## 1st Qu. 1.9999571 3.9999115
## Median 2.0000000 4.0000000
## Mean 2.0317905 3.9606228
## 3rd Qu. 2.0000751 4.0000418
## Max. 4.2001812 4.9787222
## N 201.0000000 201.0000000
## sd 0.3138139 0.3489079## Default plot method
plot(out)
## User specified plotting
matplot(out[ , 1], out[ , 2:3], type = "l", xlab = "time", ylab = "Conc",
main = "Lotka-Volterra", lwd = 2)
legend("topright", c("prey", "predator"), col = 1:2, lty = 1:2)
Example: Substrate Producer/Consumer
SPCmod <- function(t, x, parms, input) {
with(as.list(c(parms, x)), {
import <- input(t)
dS <- import - b*S*P + g*C # substrate
dP <- c*S*P - d*C*P # producer
dC <- e*P*C - f*C # consumer
res <- c(dS, dP, dC)
list(res)
})
}
## The parameters
parms <- c(b = 0.001, c = 0.1, d = 0.1, e = 0.1, f = 0.1, g = 0.0)
## vector of timesteps
times <- seq(0, 200, length = 101)
## external signal with rectangle impulse
signal <- data.frame(times = times,
import = rep(0, length(times)))
signal$import[signal$times >= 10 & signal$times <= 11] <- 0.2
sigimp <- approxfun(signal$times, signal$import, rule = 2)
## Start values for steady state
xstart <- c(S = 1, P = 1, C = 1)
## Solve model
out <- ode(y = xstart, times = times,
func = SPCmod, parms = parms, input = sigimp)
## Default plot method
plot(out)
## User specified plotting
mf <- par(mfrow = c(1, 2))
matplot(out[,1], out[,2:4], type = "l", xlab = "time", ylab = "state")
legend("topright", col = 1:3, lty = 1:3, legend = c("S", "P", "C"))
plot(out[,"P"], out[,"C"], type = "l", lwd = 2, xlab = "producer",
ylab = "consumer")
par(mfrow = mf)Example: SIR Modeling
A simple example is provided by https://www.r-bloggers.com/2017/01/sir-model-with-desolve-ggplot2/
# definition of the differential equation
sir_diffeq <- function(Time, State, Pars) {
with(as.list(c(State, Pars)), {
dSusceptible <- -beta * Infected * Susceptible
dInfected <- beta * Infected * Susceptible - gamma * Infected
dRecovered <- gamma * Infected
return(list(c(dSusceptible, dInfected, dRecovered)))
})
}
# solver for differential equation
sir_ode <- function(beta = 1, gamma = .2,
initial_state = c(Susceptible = 1-1e-6, # initial state defined
Infected = 1e-6,
Recovered = 0),
times = seq(0, 100, by = 1)) {
pars <- c(beta = beta, gamma = gamma)
# check that initial state has the correctly named variables
ode(initial_state, times, sir_diffeq, pars)
}# Parameter values
beta <- 1
gamma <- .2
out <- sir_ode(beta = beta, gamma = gamma)
out %>% as.data.frame() %>%
pivot_longer(Susceptible:Recovered, names_to = "state", values_to = "prop") %>%
ggplot(aes(time, prop)) +
geom_line(aes(color = state)) +
theme_minimal() +
labs(title = sprintf("SIR, beta = %g, gamma = %g", beta, gamma))
Example: Rossler attractor
This example is a non-linear ordinary differential equation system with chaotic dynamics
\[ \begin{aligned} \frac{dx}{dt} &= -y - z \\ \frac{dy}{dt} &= x + ay \\ \frac{dz}{dt} &= b + z(x -c) \\ \end{aligned} \]
equations x and y are linear which gives the attractor some regularity
rosslerode <- function(t, state, parms) {
with(as.list(state), {
dx <- -y - z
dy <- x + .2 * y
dz <- .2 + z * (x - 5)
return(list(c(dx, dy, dz)))
})
}
yini <- c(x=1, y = 0, z = .9)
times <- seq(0, 100, .05)
out <- ode(times = times, y = yini, func = rosslerode, parms = NULL)
scatterplot3d::scatterplot3d(out[,"x"], out[,"y"], out[,"z"], pch = 19, cex.symbols = .5)
rossler_plotly <- plot_ly(out %>% as.data.frame(), x = ~x, y=~y, z= ~z,
marker = list(color= ~time, colorscale = c("#FFFFFF", "#683531")))
rossler_plotly## No trace type specified:
## Based on info supplied, a 'scatter3d' trace seems appropriate.
## Read more about this trace type -> https://plotly.com/r/reference/#scatter3d## No scatter3d mode specifed:
## Setting the mode to markers
## Read more about this attribute -> https://plotly.com/r/reference/#scatter-modeExample: Van der Pol Equation
This is a good example of a system that can be either stiff or non-stiff.
- For large values of \(\mu\), it’s stiff.
- For small values of \(\mu\), it’s non stiff
The 2nd order ODE can be transformed into a 2 first order equations:
\[ \begin{aligned} y'' - \mu(1 - y^2)y' + y &= 0 \\\\ x' &= y\\ y' &= \mu(1 - x^2)y - x \end{aligned} \]
vdpol <- function(t, state, mu) {
with(as.list(state), {
list(c(y, mu * (1 - x^2) *y - x))
})
}
yini <- c(x = 2, y = 0)
vdpol_stiff <- ode(yini, times = 0:3000, func = vdpol, parms = 1000)
vdpol_nonstiff <- ode(yini, times = seq(0, 30, .01), func = vdpol, parms = 1)
plot(vdpol_stiff, type = "l", main = "IVP ODE, stiff", xlim = c(0, 3000), which = "x")
plot( vdpol_nonstiff,type = "l", main = "IVP ODE, nonstiff", xlim = c(0, 30), which = "x")